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Abstract 

A characteristic feature of functional data is the presence of time variability in addition 
to amplitude variability. The existing functional regression methods do not handle time 
variability in an explicit and efficient way. In this paper we introduce a functional regres- 
sion method that incorporates time warping as an intrinsic part of the model. The method 
achieves good predictive power in a parsimonious way, and allows for unified statistical 
inference of time and amplitude variability. The properties of the estimators are studied 
by simulation, and an application to the modeling of ground-level ozone trajectories is 
presented. 

Key Words: Functional Data Analysis; Random-Effect Models; Registration; Spline 
Smoothing; Time Warping. 



1 Introduction 

Many statistical applications involve modeling curves as functions of other curves. For 
example, the evolution of CD4 cell counts of HIV patients can be modeled as a function of 
viral load (Liang et al. 2003, Wu and Liang 2004, Wu and Miiller 2010); gene expression 
profiles of an insect at the pupa stage can be modeled as functions of gene expression pro- 
files at the embryonic stage (Miiller et al. 2008); the trajectory of systolic blood pressure 
over the years can be explained to some extent by the trajectory of body mass index (Yao 
et al. 2005). 

Functional regression, as used in the above papers, is a rather straightforward exten- 
sion of multivariate linear regression to functional data (Ramsay and Silverman 2005, 
ch. 16). The method consists of finding the principal components of the explanatory and 
the response functions, and fitting a multivariate linear model on the respective component 
scores. Different procedures may be used to estimate the principal components (depending 
e.g. on whether the sample curves are pre-smoothed or the raw data is used), but always 
the main idea is to explain the modes of variation of the response curves in terms of the 
modes of variation of the explanatory curves. This makes the functional regression model 
easy to interpret if the principal components are themselves easy to interpret, but this is 
not always the case. 

Recent developments in functional linear regression have focused on theoretical as- 
pects such as convergence rates (Cai and Hall 2006, Hall and Horowitz 2007, Crambes et 
al. 2009), on sparse longitudinal data (Yao et al. 2005), and on the interpretability of the 
estimators (James et al. 2009). But a distinct feature of functional data that has not yet 
been properly investigates in the regression context is time variability. It is well known 
that functional samples normally present "horizontal" variability of the features (e.g. lo- 
cation of the peaks) in addition to "vertical" variability (e.g. amplitude of the peaks), and 
that a preliminary alignment (or warping) of the curves can considerably reduce the num- 
ber of principal components that are necessary to explain amplitude variability, and at 
the same time, make the principal components easier to interpret (Ramsay and Silver- 
man 2005, ch. 7). Many warping methods have been proposed in recent years (Gervini 
and Gasser (2004, 2005), James (2007), Kneip et al. (2000), Kneip and Ramsay (2008), 
Liu and Miiller (2004), Ramsay and Li (1998), Tang and Miiller (2008, 2009), Wang and 



Figure 1 : Ozone Example. Daily trajectories of ground-level concentrations of (a) oxides 
of nitrogen and (b) ozone in the city of Sacramento in the Summer of 2005. 

Gasser (1999)). But classical functional linear regression neither explicitly models time 
variability nor eliminates it by warping. Since time variability tends to spread the features 
of the curves over wide time ranges, a large number of principal components are typically 
needed to provide a good fit to the data, and the corresponding regression model is often 
overparametrized. 

Synchronizing the curves before fitting the regression model would make it more par- 
simonious and more interpretable. But once time variability is removed, the potentially 
useful information contained in the warping process would be lost. Moreover, since the 
model would ignore time variability in the response curves, its predictive power would be 
low. One possible solution would be to use two regression models, one for the amplitude 
components and one for the warping components, but statistical inference from such a 
two-stage approach would be cumbersome, if at all possible. 

In this paper we introduce a functional regression model that incorporates time warping 
explicitly into the model. This is accomplished via a random-effects simultaneous linear 
model on the warping parameters and the principal component scores. This approach 
unifies, but does not confound, the treatment of amplitude and time variability, resulting 
in a regression model that is more efficient, parsimonious and interpretable than classical 
functional linear regression. To better illustrate the potential advantages of incorporating 
time variability into the model, consider the sample curves in Figure [Q which will be 



analyzed in more detail in Section |5J The curves in Figure HJa) are daily trajectories of 
oxides of nitrogen in the city of Sacramento, California, during 52 summer days of 2005. 
Figure HJb) shows the corresponding trajectories of ozone concentration. The goal is to 
predict ozone concentration trajectories from oxides-of- nitrogen trajectories. Both types 
of curves follow regular patterns; in particular, they show peaks that vary not only in 
amplitude but also in timing. It is reasonable to hypothesize not only that a large peak in 
oxides of nitrogen will be followed by a large peak in ozone concentration, but also that 
an early peak in oxides of nitrogen will be followed by an early peak in ozone level, and 
vice- versa. Therefore, both amplitude and timing of the features of the explanatory curves 
can provide useful information for predicting the response curves. This, in fact, will be 
shown in Section |5]for these data. 

The paper is organized as follows. The warped regression mode is presented in Section 
[2] Estimation and prediction is discussed in Section [3j The properties of the estimator are 
studied by simulation in Section HI and the application of the new method to the ozone 
data is presented in Section |5] A technical supplement and Matlab/Fortran programs are 
available on the author's website. 

2 The Warped Functional Regression Model 

Consider a sample of pairs of functions (xi, yi), ..., (x n , y n ), where Xi(s) is the covariate 
and yi(t) the response. We assume Xi : 5? — > K and yi : 3? — > K with 5? and ST closed 
intervals in R. For simplicity we will refer to s and t as "time", although they could be any 
one-dimensional variable (such as wavelength, if the data are spectral curves.) For most 
applications 5? will be equal to 2f , but this is not strictly necessary and it is not assumed. 
Normally the sample functions Xi(s) and yi(t) are observed at discrete time points 
{ s ij '■ J — I; • • • j v u} an d {Uj : j = 1, . . . , z/ 2 j}, respectively. In many cases the grids 
are dense enough as to allow individual smoothing of the curves, but in other cases the 
number of observations per curve is low and the time points are sparse. This is a more 
typical situation in longitudinal data analysis, and it is the kind of data we will focus on 
in this paper. The methods developed here are impractical for u u or u 2 i greater than 30 or 
so, although they can still be applied to smooth curves if the curves are down-sampled and 
some noise is added. 



Let us assume then that the raw data are vectors (xi, yi), . . . , (x n , y n ) such that, for 
each subject i = l,...,n, 

Xij = Xi(Sij) + €ij, €ij ~ JV(0, <J 2 £ ), j = 1, . . . , U a , 

Uij = ViiUj) + Vij, Vij ~ N(0, crj), j = 1, . . • , z/ i2 . 
The curves {xj(s)} and {^(t)} can be thought of as compound processes 



Xi{s) = x*(ur\ s )), (1) 

Vi(t) = tfCCT 1 ®), ( 2 ) 



where {x*(s)} and {?/*(£)} present only amplitude variability in the sense described below, 
and {ui(s)} and {d(t)} are warping functions (i.e. U{ : 5? — >■ J^ and (^ : 5^ — y ^ are 
strictly monotone increasing) accounting for the time variability. Of course, for any given 
UJi(s) and d(t) one can always make (Q} and © hold just by defining s*(s) = Xi(ui(s)) 
and y*(t) = yi(d(t)), but the goal is to find families of functions {ui(s)}, {d(t)}, {x*(s)} 
and {y*(t)} that live in low-dimensional spaces. We assume then that {x*(s)} and {y*(t)} 
follow Karhunen-Loeve-type decompositions 



pi 



x*(s) = /j, x (s) + ^2u ik (f) k (s), u ik ~ JV(0, A fc ) independent, (3) 

fc=i 

y*(t) = fi y (t) + ^2v il i/) l (t), vu ~ JV(0, 7j) independent, (4) 

z=i 

where {0 fe (s)} and {^(t)} are orthonormal basis functions in their respective spaces 
L 2 (,y) and L 2 (^), and the dimensions p\ and p 2 are low. For identifiability purposes 
we assume Ai > • • ■ > A Pl > and j x > ■ ■ ■ > 7 P2 > 0. The basis functions, which we 
will call KL-components for short, will be estimated from the data (see Section[3]) 

For the warping functions we will assume semiparametric models based on interpo- 
lating Hermite splines. Although other models are possible, using e.g. integrated splines 
(Ramsay, 1988) or monotonically-constrained B-splines (Brumback and Lindstrom, 2004), 
Hermite splines offer flexibility, parsimony and interpret ability of the parameters at the 
same time. Given r knots (or landmarks) in y = [a,b], a < t 01 < ■ ■ ■ < r 0r < b, and 



correspondingly r individual landmarks a < rn < ■ ■ ■ < r ir < b for each i = 1, . . . , n, it 
is always possible to find a strictly increasing piecewise cubic polynomial Ui(s) such that 
Wt(roj) = Tij for all j. Moreover, we can express Ui as 

r+l r+1 

3=0 3=0 

where {aj(s; t )} and {/3,-(s; r )} are polynomial functions that only depend on t , and 
the dijS (which satisfy dy = ^(r^)) can be deterministically obtained from the r^-s via 
the Fritsch-Carlson algorithm (Fritsch and Carlson, 1980). (Details are given in the tech- 
nical supplement.) The individual landmarks r^ will be treated as random effects. The 
knots r are chosen by the user. For instance, if the sample curves show a prominent peak 
(as in Figure [Q), one can take r as the (approximate) average peak location, and the cor- 
responding Ti would be interpreted as the peak location of the zth curve (the fundamental 
difference with traditional "landmark registration" is that we do not estimate the individual 
TjS but rather treat them as random effects.) Since the number of salient features of the 
curves is typically small, the dimension r will also be small. 

Due to the monotonicity restriction on the elements of Tj, it is not easy to specify 
a simple distribution for them. So we will proceed as in e.g. Brumback and Lindstrom 
(2004) and use the Jupp transform (Jupp 1978): 

g0=bgf r *+ 1 ~ T * Y J = l,...,r, 
\ T ij r i,j-i/ 

where r, = a and r, ; r+1 = b. The vector t is similarly transformed into 6 . The 6s 
are now unrestricted vectors in IR r , so we can reasonably assume a Normal distribution for 
them. To summarize, our warping model can be expressed as: 

wr 1 ^) = 9x(s;9 X i,6 x0 ), xi ~ N(0 xO ,'E ex ), (6) 

Ci\t) = g y (t;0 yt ,0 y0 ), yi ~N{e y0 ,Ve y ), (7) 

for fixed functions g x and g y whose exact form does not matter. The covariance matrices 
Tt 9x and S 0y will be estimated from the data (see Section |3}. 

Our proposed regression model assumes the response effects {v,} and {6 yi } are linear 



functions of the covariate effects {uj} and {0 xi }. Specifically, let 
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where A = diag(Ai, . . . , A Pl ) and T = diag(7 1; . . . , 7 ). We assume 



z, = /x 2 + A (wi - /xj + e i7 e< ~ JV(0, E e ), 



(8) 



with S e diagonal and A a (p 2 + ^2) x (pi + r i) matrix, which will be the main parameter 
of interest. 

For interpretability it is convenient to split A into four blocks corresponding to the 
subvectors that make up Wj and z$: 



A 



An A12 

A21 A 2 2 



Since T, z = AS W A T + S e , the condition that T be diagonal translates into the restriction 
that [An, AijEtJAii, Ai 2 ] T be diagonal. 



3 Estimation 

The means and KL-components in ([3]) and Q are estimated from the data using B-splines. 

Letb^s) = (6^1(5),..., fe xgi (s)) T beaB-spline basis inL 2 (^)andbj / (t) = (b yl (t), . . . ,b yq2 (t)) T 
a B-spline basis in L 2 (3?). We assume jjl x (s) = b^(s)m x , [i y (t) = hy(t)niy, <f> k (s) = 



h?(s)c k aadipi(t) = b y (t)di for certain parameters m x , m y , {c k } and {d^}. Let C = [c v . . . , c Pl ] 
and D = [d 1; . . . , d p2 \. Then the orthogonality restrictions on the KL-basis translate into 
the restrictions C T J X C = I P1 and D T J J/ D = I P2 , where J x = Jb x (s)b^(s)ds and 
J y = Jb y (t)bl(t)dt. 

To summarize, our model has the following parameters: A, S e , S„, m^, m^, C, 
D, <Tg and crjL In addition, the following constants must be chosen by the user: the model 
dimensions pi, p 2 , r\ and r 2 ; the warping knots t x0 and T y0 , which determine 6 x0 and y0 ; 
and the spline basis parameters (number and placement of the knots, and spline order.) In 
our experience, cubic splines with equally-spaced knots work well for practical purposes; 
the number of knots can be adjusted depending on the smoothness of the KL-components. 
The model dimensions p±, p 2 , r\ and r 2 could be chosen by Akaike-type criteria, but more 
realistically, they will be suggested by the number of features in the sample curves and 
by the marginal tests of significance of the coefficients of A, as it is customarily done in 
linear regression. 

The parameters can be estimated by maximum likelihood using the EM-algorithm. 
The complete-data likelihood is derived as follows. Let Sj = (sn, . . . , Si Uil ) T and tj = 
(tn, . . . ,tiu i2 ) T be the individual time grids in vector form, s* = g x (si] 6 xi , xO ) and 
t* = g y (ti, Byi, y o) the corresponding warped grids, and B* xi = [6 x i(s*), . . . , b xqi (s*)] 
and B*j = [b y i(t*), . . . , b yq2 (t*)] the corresponding warped basis functions (note that s* 
and t*, and thus B*,- and B* i5 are functions of xi and 6 yi and therefore random quantities.) 
Then we have 

Xi|wj ~ ^(B^m^ + B^Cu^cr^J, 

yi\zi ~ ^(B^m^ + B^Dv^aJl^), 

Zi|wi ~ N(n z + A(wi-/*J,S e ), 

Wj ~ iV(/* w ,S w ). 

Treating the random effects (wj, Zj) as missing data, the complete-data density factorizes 

as 

f(xi,yi,Wi,Zi) = f(x i ,y i \w i ,z i )f(zi\w i )f(wi) 

= f(^i\w i )f(y i \z i )f(z i \w i )f(w i ), 



from which updating equations for the parameters can be derived. The restrictions on C, 
D, and A do pose some difficulties, so the implementation of the EM algorithm is not 
trivial; see the technical supplement for details. Our implementation estimates the mod- 
els sequentially, in increasing order of complexity: we start with a warping-only model 
(i.e. pi = p 2 = 0), and add KL-components one at the time, using the estimators of the 
previous model as initial estimators of the larger model, until the desired p\ and p-2 are 
reached. Matlab/Fortran programs implementing this algorithm are freely available on the 
author's website. 

In addition to parameter estimation, prediction of a new response given a new covariate 
is often of interest. This is straightforward for our model. Interestingly, since no assump- 
tions on the time grids are made we can really predict "the future", that is, we can predict 
the response curve at later time points than those available for the explanatory curve. Let 
x be the covariate data observed at a time grid s . The unobserved random effect w is 
estimated by w = ^(wqIxq), where E denotes expectation computed with the estimated 
parameters. More explicitly, 

- _ Jw/(x |w)/(w)dw 

Wg — , 

f/(x |w)/(w)dw 

where the integrals are approximated by Monte Carlo integration. The predictor of z is 

z = -E(z |x ) = (i z + A (w — A™)- Once 2q = [vq , 6 y0 \ is obtained, the predicted am- 
plitude and warping components of the response curve are i)o(t) = fi y {t) + Y7i=i voi'<fii(t) 
and £o (t) = g y (t; 6 y0 , y0 ), respectively, from which y (t) = yg(C (*)) is obtained. 



4 Simulations 

We ran a number of simulations to study two aspects of the proposed method: estimation 
accuracy (as a function of the sample size, the grid size and the spline-basis dimension), 
and prediction accuracy (particularly compared with common functional linear regres- 
sion). 

To study estimation accuracy, we simulated a simple model with one component and 
one warping knot. Specifically, if (p(s,fi,a) denotes the N(fi,a 2 ) density function, we 
took fi x (s) = .6y?(s, .3, A)+A(p(s, .6, .1), ^(s) = cp(s, .3, .1) / '1.6796, n y (t) = .6y?(£, .5, .1) + 

8 



Figure 2: Simulation Results. A sample of ten simulated curves: (a) explanatory curves 
and (b) the corresponding response curves. 

Aip(t, .8, .1) and i> x (t) = <p(t, .5, .1)/1.6796, with s and t in [0, 1]. The warping knots 
were r x0 = .3 and r y0 = .5. Thus, although the mean functions have two peaks, hor- 
izontal and vertical variability is concentrated on the main peak. The regression matrix 
A was taken as the identity matrix, and the other parameters were £„, = diag(.2 2 , .2 2 ), 
S e = diag(.05 2 , .05 2 ), and a e = a r] = .1. The time grids s$ and tj were equally- spaced 
points in [0, 1]. A sample of 10 curves is shown in Figure [T]for illustration. 

For estimation of the functional parameters we used cubic B-splines with equally 
spaced knots; we considered two scenarios, which represent mild oversmoothing and mild 
undersmoothing: 5 knots (resulting in 9 basis functions), and 10 knots (14 basis functions). 
Two sample sizes were considered, n = 50 and n = 100, and two grid sizes, m — 15 and 
m = 30. Each of these eight combinations were replicated 500 times. For a functional 
parameter (/, say), we defined as "bias" the quantity {f(E(f) — f) 2 } 1 ^ 2 and as "standard 
deviation" the quantity {J V(f)} 1 ^ 2 (for scalars, the usual definitions were used.) 

The results are summarized in Tables [TJand |2] We see that the most important factor 
that determines the performance of the estimators is the number of knots; which is to say, 
the bias of the estimators of the functional parameters. One might have thought that using 
too many knots for a small sample size would have led to a large estimation error due to the 
larger variance, but that was not the case: the bias clearly dominates, and the significant 
bias reduction attained by increasing the number of knots compensates for the variance 



Parameter 

a n 

0-12 
021 

«22 
//,. 



/;/ 



n = 

= 15 

.13 

.34 

.21 

4.01 

2.23 

.64 

2.18 

.56 



50 
m 



knots 

= 30 

.15 

.38 

.08 
4.52 
2.34 

.58 
2.29 

.59 



in 



n = 

= 15 

.13 

.38 

.21 

4.09 

2.23 

.65 

2.25 

.56 



100 

m ■ 



= 30 

.14 

.34 

.01 

4.46 

2.34 

.58 

2.39 

.58 



n = 50 
m — 15 m 
.00 
.01 
.09 
.13 
.09 
.09 
.73 
.80 



knots = 10 

= 30 m 
.09 
.07 
.21 
.17 
.10 
.09 
.81 
.74 



n = 100 
= 15 m = 30 



.05 
.02 
.13 
.15 
.08 
.08 
.79 
.78 



.07 
.00 
.19 

.17 
.09 
.08 

.84 

.77 



Table 1: Simulation Results. Bias (x 10) of the estimators. 



Parameter 

a n 

Ol2 
021 

«22 

»x 

U y 






n = 50 
m — 15 m 
.79 
.74 
1.28 
.80 
.43 
.42 
.91 
.78 



knots = 

= 30 m 

.78 

.68 
1.28 

.65 

.40 

.40 

.81 

.76 



n = 100 
= 15 m = 30 



.57 
.53 
.91 
.60 
.31 
.28 
.66 
.59 



.50 
.46 
.84 
.46 
.29 
.27 
.61 
.49 



n = 50 
m — 15 m 

.58 

.58 

.99 

.56 

.46 

.55 
1.17 
1.25 



knots = 10 

= 30 m 
.58 
.56 
.94 

.54 
.47 
.55 
.92 
1.12 



n = 100 
= 15 m = 30 



.42 
.40 
.59 
.40 

.33 
.38 

.74 
.84 



.37 
.35 
.63 
.37 
.32 
.37 
.64 
.71 



Table 2: Simulation Results. Standard deviation (x 10) of the estimators. 



increase. Since the model "borrows strength" across curves to estimate the structural pa- 
rameters, it is possible to estimate a total of 52 spline coefficients with only 50 curves in a 
reasonably accurate way. The reason the estimation bias of the functional parameters has 
such a large effect on A seems to be attributable to the fact that more accurate estimation 
of the peaks translates into a better alignment of the curves, and consequently a better 
estimation of the regression coefficients. This can be gleaned from individual plots of the 
estimators, which are not shown here for reasons of space. The grid size m does not have 
much influence on the accuracy of the estimators. 

To assess prediction accuracy we simulated data from the above model with n = 50, 
m = 20, and four different combinations of warping-effect variance H$ x and random- 



10 





S = 

y 


.0425 


S e = 

y 


.1625 




r 2 = .94 


rp£ 


= .70 


r 2 = .98 


r 2 = .70 


FLR-1 


.277 




.299 


.322 


.373 


FLR-4 


.175 




.211 


.247 


.334 


FLR-9 


.156 




.204 


.194 


.311 


WFR 


.143 




.203 


.145 


.317 



Table 3: Simulation Results. Prediction errors of ordinary functional linear regression 
(FLR) of various dimensions (1, 4 and 9) and warped functional regression (WFR). 

error variance £ e; 22- As before, we considered Y>e x = -2 2 and £ e ,22 = -05 2 , and also a 
larger variance T,g x = -4 2 with £ ej2 2 = -05 2 ; the corresponding y-warping variances are 
S e = .0425 and S e = .1625, and the corresponding coefficients of determination r 2 = 
E^/Efl are .94 and .98. Since these r 2 s are rather high, we also simulated two additional 
models with E 6y = .0425 and Y> dy = .1625 but with r 2 = .70 (for which Y> dx = .17 2 and 
E 6) 22 = -ll 2 , and Eg x = .34 2 and E e2 2 = -22 2 , respectively.) A total of 500 replications 
for each model were run. For each replication, a further test set of n* = 50 observations 
was generated and the responses predicted as explained in Section [3l prediction error 
was measured by the root average squared error, {$^"=1 E(\\yi ~ 9i\\ )/rnn*} 1 / 2 . For 
comparison we also fitted a classical functional linear regression, as in e.g. Miiller, Chiou 
and Leng (2008), with the difference that the KL-components were estimated via random- 
effect models (James, Hastie and Sugar, 2000) rather than by kernel smoothing. Since the 
number of KL-components determines the classical functional regression estimator, we 
considered three cases: 1, 2 and 3 components for both for the x's and the y's, so that the 
model dimensions were 1, 4 and 9, respectively. 

Prediction errors are reported in Table For classical functional linear regression, the 
error decreases as the model dimension increases. This is not unexpected, since horizontal 
variability can always be modeled as vertical variability provided a large number of com- 
ponents are used. The real question is how overparameterized ordinary linear regression 
needs to be in order to attain a prediction error comparable to a low-dimensional warped 
regression. As Table |3] shows, for large r 2 not even a 9-dimensional classical regression 
model attains an error comparable to the 4-dimensional warped regression model. For 
r 2 = .70 the 9-dimensional model does attain a prediction error comparable to warped 
regression, but from a qualitative point of view warped regression provides more neatly 
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interpretable amplitude components and warping landmarks, whereas ordinary linear re- 
gression provides KL-components that are hybrids of amplitude and time variability, to 
the point that the third component is almost impossible to interpret. 

5 Application: Modeling Ground-Level Ozone Concen- 
tration 

Ground-level ozone is an air pollutant that is known to cause a number of serious health 
problems. Unlike other pollutants, ozone is not emitted directly into the air but is formed 
as a result of complex chemical reactions, including among other factors, volatile organic 
compounds and oxides of nitrogen in the presence of heat and sunlight. Volatile organic 
compounds and oxides of nitrogen are emitted from a variety of sources, such as mo- 
tor vehicles and combustion engines, chemical plants, power plants, and other industrial 
sources. Modeling ozone formation has been an active topic of air-quality studies for many 
years, and the constant tracking of air-pollution levels in most cities of the United States 
has provided researchers with enormous databases. 

In this article we will make use of the California Environmental Protection Agency 



database, available at http://www.arb.ca.gov/aqd/aqdcd/aqdcddld.htm. Hourly concentra- 



tions of many pollutants at many locations in California are available for the years 1980- 
2009. For the purposes of this section, which is just to give an example of application of 
warped functional regression, we will focus on the relationship between the trajectories 
of oxides of nitrogen (NOx) and ozone (03) in the city of Sacramento (site 3011 in the 
database), in the Summer of 2005. We skipped weekends and holidays because the NOx 
and 03 levels are substantially lower than for weekdays, due to the reduced car traffic and 
industrial activity. There were also a number of outlying trajectories that were eliminated, 
which left a sample of 52 days between June 6 and August 26. These are shown in Figure 

Both NOx and 03 trajectories show simple regular patterns. NOx curves tend to peak 
around 8am and 03 around 3pm. Therefore, we fitted warped regression models with 
one warping knot. We tried several combinations of t x0 and r y0 around 8am and 3pm, 
respectively, and in all cases the estimated means turned out to have peaks at 7am for 



12 



NOx and 2pm for 03 (which shows that warped regression is robust to specification of 
r x0 and r y0 , at least within a reasonable range); the estimators reported here correspond 
to r x0 = 7 and r y0 = 14. As spline basis functions we used cubic B-splines with 7 
equispaced knots (one every 3 hours); we also tried 10 knots but the results were not 
substantially different. Three models were fitted: (i) one KL-component for each of x 
and y, (ii) two KL-components for x and one for y, and (Hi) one KL-component for x 
and two for y. The log-likelihood values were 44.44, 45.21 and 52.04, respectively. The 
second model did not seem to represent much of an improvement over the first one, so we 
discarded it. For models (i) and (Hi) the estimated regression coefficients and the bootstrap 
standard deviations (based on 200 resamples) were: 



A 



.726 .085 

.192 .439 

.355 .120 

.010 .015 

.180 .540 



st.dev.(A) 



st.dev.(A) 



.070 .021 

.079 .056 

.081 .065 

.036 .097 

.057 .114 



While all regression coefficients of model (i) are statistically significant, for model (ii) 
none of the coefficients of the second y-component are significant, so we chose model (i). 
The other parameter estimators for this model are: 



.0206 -.0330 
-.0330 .0863 



.0021 
.0042 



.0094 
-.0050 



-.0050 
.0160 



a 



.351 and a v = .021. To help interpret the KL-components, Figure |3£a) shows \i x 
and fj, x ± 2^/\[(j) 1 and Figure |3{b) shows // and [i y ± l^fx^x- Both are shape compo- 
nents with basically the same interpretation: curves with positive scores will tend to have 
sharper features than the mean while curves with negative scores will tend to be flatter 
than the mean. The diagonal coefficients of A being positive indicates that the component 
scores Ui and t)j are positively correlated, as Figure |3jc) shows, and the warping landmarks 
f xi and T y i, which for this model can be roughly interpreted as the peak locations, are also 
positively correlated, as Figure |3jd) shows. Note that the warping effects and the compo- 
nent scores are highly correlated: from the above estimates we have coYr(9 xi , u») = —.78 
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Figure 3: Ozone Example. Warped Functional Regression fit. (a) Log-NOx mean (solid 
line), and mean plus (dashed line) and minus (dotted line) two standard deviation times 
the KL component; (b) same as (a) for the square root of 03; (c) NOx versus 03 KL- 
component scores; (d) NOx versus 03 estimated landmarks. 
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Figure 4: Ozone Example. Weekday ozone trajectories in July 2006 [(a)], and predicted 
trajectories using warped functional regression [(b)], 4-dimensional common functional 
regression [(c)], and 9-dimensional common functional regression ([d]). 

and corr(9 yi , Vi) = —.40 (note that the Jupp transformation is a decreasing function of the 
knots, so the correlations between f xi and Ui and between r yi and Vi are actually posi- 
tive). Therefore, curves with late peaks tend to have sharper features than the mean, and 
vice-versa. 

It is also of interest to investigate to which extent this model can be used to predict 03 
curves in subsequent years. So we chose a test sample of NOx curves for the month of 
July of 2006, a total of 19 after eliminating some days with missing data. For compari- 
son, we also computed the predictors from ordinary functional linear regression with 1, 2 
and 3 KL-components (as in Section (U). The resulting root mean prediction errors were 
comparable: .030 for the warped linear regression and ordinary functional regression with 
1 and 2 components, and .029 for functional regression with 3 components. The corre- 
sponding curves are shown in Figure |4] Perhaps the difference between both methods is 
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not very pronounced because there is not much time variability in the response peaks, but 
from a qualitative point of view, we see that warped regression provides better prediction 
of the dynamics of the peak than the 4-dimensional ordinary linear regression; only the 
9-dimensional ordinary linear regression model is able to predict dynamic variation of the 
peaks as well as the warped regression model. This is in agreement with the simulation 
results of Section 0] 

Acknowledgement 

This research was partially supported by NSF grant DMS 10-06281. 

References 

Brumback, L.C. and Lindstrom, M.J. (2004). Self modeling with flexible, random time 
transformations. Biometrics 60 461-470. 

Cai, T and Hall, P. (2006). Prediction in functional linear regression. The Annals of 
Statistics 34 2159-2179. 

Crambes, C, Kneip, A., and Sarda, P. (2009). Smoothing splines estimators for func- 
tional linear regression. The Annals of Statistics 37 35-72. 

Fritsch, F.N. and Carlson, R.E. (1980). Monotone piecewise cubic interpolation. SIAM 
Journal of Numerical Analysis 17 238-246. 

Gervini, D. and Gasser, T (2004). Self-modeling warping functions. Journal of the Royal 
Statistical Society (Series B) 66 959-971. 

Gervini, D. and Gasser, T (2005). Nonparametric maximum likelihood estimation of the 
structural mean of a sample of curves. Biometrika 92 801-820. 

Hall, P. and Horowitz, J. L. (2007). Methodology and convergence rates for functional 
linear regression. The Annals of Statistics 35 70-91 . 



16 



James, G.M. (2007). Curve alignment by moments. The Annals of Applied Statistics 1 
480-501. 

James, G., Hastie, T. G. and Sugar, C. A. (2000). Principal component models for sparse 
functional data. Biometrika 87 587-602. 

James, G., Wang, J. and Zhu, J. (2009). Functional linear regression that's interpretable. 
The Annals of Statistics 37 2083-2108. 

Jupp, D. L. B. (1978). Approximation to data by splines with free knots. SI AM J. Numer. 
Anal. 15 328-343. 

Kneip, A., Li, X., MacGibbon, B. and Ramsay, J.O. (2000). Curve registration by local 
regression. Canadian Journal of Statistics 28 19-30. 

Kneip, A. and Ramsay, J.O. (2008). Combining registration and fitting for functional 
models. Journal of the American Statistical Association 103 1 155-1 165. 

Liang, H., Wu, H., and Carroll, R. J. (2003). The relationship between virologic and im- 
munologic responses in AIDS clinical research using mixed-effects varying-coefficient 
models with measurement error. Biostatistics 4 297-312. 

Liu, X. and Miiller, H.-G. (2004). Functional convex averaging and synchronization for 
time- warped random curves. Journal of the American Statistical Association 99 
687-699. 

Miiller, H.-G., Chiou, J.-M., and Leng, X. (2008). Inferring gene expression dynamics 
via functional regression analysis. BMC Bioinformatics 9 60. 

Ramsay, J.O. (1988). Monotone regression splines in action (with discussion). Statistical 
Science 3 425-461. 

Ramsay, J.O. and Li, X. (1998). Curve registration. Journal of the Royal Statistical 
Society (Series B) 60 351-363. 

Ramsay, J.O. and Silverman, B. (2005). Functional Data Analysis (Second Edition). 
Springer, New York. 



17 



Tang, R., and Miiller, H.-G. (2008). Pairwise curve synchronization for functional data. 
Biometrika 95 875-889. 

Tang, R., and Miiller, H.-G. (2009). Time-syncronized clustering of gene expression 
trajectories. Biostatistics 10 32-45. 

Yao, R, Miiller, H.-G. and Wang, J.-L. (2005). Functional linear regression analysis for 
longitudinal data. The Annals of Statistics 33 2873-2903. 

Wang, K. and Gasser, T (1999). Synchronizing sample curves nonparametrically. The 
Annals of Statistics 27 439-460. 

Wu, H. and Liang, H. (2004). Backfitting random varying-coefficient models with time- 
dependent smoothing covariates. Scandinavian Journal of Statistics 31 3-19. 

Wu, S. and Miiller, H.-G. (2010). Response-additive regression for longitudinal data. 
Unpublished manuscript. Available at http://anson.ucdavis.edu/~mueller/rare8.pdf. 



18 



(a) 



0.3 










0.3 


0.25 










0.25 


co 0.2 
O 










0.2 


ct 0.15 

C/3 










0.15 


0.1 










0.1 


0.05 


P^f^ 








0.05 


( 


) 5 


10 


15 


20 






t (hour) 












(c) 








0.3 










0.3 


0.25 










0.25 


0.2 










0.2 


0.15 










0.15 


0.1 










0.1 


0.05 










0.05 




10 



15 



20 





10 15 

t (hour) 




10 15 

t (hour) 



20 





0.2 0.4 0.6 0.8 




10 15 

t (hour) 















(b) 










0.3 












/ 

/ 
I 


*-■ 


\ 


\ 


\ 


0.25 




















\ 
\ 
























8 °- 2 












I / 








\ \ 














l / 




























w 0.15 










/l 


1/ . ■ ' ' 








\\ 
. \\ 


0.1 


\ 


\ 


^ 


/ 


/ 












0.05 




















- 



5 10 15 

t (hour) 



20 















(c) 






0.3 






























0.2 












o o 
o ° 


CD 
O 


0.1 










oo 


tf^ 


O 
CO 

I 













^°°° %s 


o 


CO 

O 


-0.1 
-0.2 


O 


o 


o 
o 


o 


% o o 


- 



(d) 



-0.4 -0.2 0.2 

NOx KL-score 



0.4 



16 



15 

ST 

<D 

n 14 
O 



13 



12 





o 


o 


" 


o o£$o 


" 


o6Hgr 


o^° 
oo° o 

3? 


- 


6to 







10 11 



x (NOx peak) 



